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DEVELOPMENT AND ANALYSIS OF A 
BLOCK-PRECONDITIONER FOR THE PHASE-FIELD CRYSTAL 

EQUATION 

SIMON PRAETORIUS+t AND AXEL VOIGT* 

Abstract. We develop a preconditioner for the linear system arising from a finite element 
discretization of the Phase Field Crystal (PFC) equation. The PFC model serves as an atomic 
description of crystalline materials on diffusive time scales and thus offers the opportunity to study 
long time behaviour of materials with atomic details. This requires adaptive time stepping and 
efficient time discretization schemes, for which we use an embedded Rosenbrock scheme. To resolve 
spatial scales of practical relevance, parallel algorithms are also required, which scale to large numbers 
of processors. The developed preconditioner provides such a tool. It is based on an approximate 
factorization of the system matrix and can be implemented efficiently. The preconditioner is analyzed 
in detail and shown to speed up the computation drastically. 
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1. Introduction. The Phase Field Crystal (PFC) model was introduced as a 
phenomenological model for solid state phenomena on an atomic scale [26, 27]. How¬ 
ever, it can also be motivated and derived through classical dynamic density functional 
theory (DDFT) [28, 65] and has been used for various applications in condensed and 
soft matter physics, see the review [31] and the references therein. Applications in¬ 
clude non-equilibrium processes in complex fluids [7, 53], dislocation dynamics [19], 
nucleation processes [12, 9, 35, 13], (dendritic) growth [32, 70, 62] and grain growth 
[ 10 ]- 

The main solution methods for the PFC model, which is a non-linear 6th order 
parabolic partial differential equation, are finite-difference discretizations and spectral 
methods, which are combined with an explicit or semi-implicit time-discretization. 
Numerical details are described in [20, 37, 38, 63, 29]. 

Recently, the PFC model has been coupled to other field variables, such as flow 
[53], orientational order [1, 54] and mesoscopic phase-field parameters [45]. This limits 
the applicability of spectral methods due to the lack of periodic boundary conditions 
in these applications. On the other hand, simulations in complex geometries have been 
considered, e.g. colloids in confinements motivated by studies of DDFT [5], crystal¬ 
lization on embedded manifolds [14, 11,4, 60] or particle-stabilized emulsions, where 
the PFC model is considered at fluid-fluid interfaces [2, 3]. These applicabilities also 
limit the use of spectral and finite-difference methods or make them sometimes even 
impossible. The finite element method provides high flexibility concerning complex 
geometries and coupling to other partial differential equations, which is the motiva¬ 
tion to develop efficient solution methods for the PFC model based on finite element 
discretizations. 

Basic steps of finite element methods include refinement and coarsening of a 
mesh, error estimation, assembling of a linear system and solving the linear system. 
Most previous finite element simulations for the PFC model [12, 14, 53] have used 
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direct solvers for the last step, which however restrict the system size due to the high 
memory requirements and only allow computations in 2D. Well-established solution 
methods for linear systems, such as iterative Krylov-subspace solvers, like CG, MIN- 
RES, GMRES, TFQMR, BiCGStab are not directly applicable for the PFC equation 
or do not converge, respectively converge very slowly, if used without or with standard 
preconditioners, like Jacobi or ILU preconditioners. 

In this paper, we propose a block-preconditioner for the discretized PFC equations 
and analyze it with respect to convergence properties of a GMRES method. We 
have organized the paper as follows: In the next section, we formulate the PFC 
model in terms of a higher order non-linear partial differential equation. Section 
3 introduces a space- and time-discretization of the model, including the treatment 
of the non-linearity. In Section 4, the preconditioner is introduced and an efficient 
preconditioning procedure is formulated. The convergence analysis of GMRES is 
introduced in Section 5 and Section 6 provides an analysis of the preconditioner in 
terms of a spectral analysis. Finally, in Section 7 we examine the preconditioner 
in numerical examples and demonstrate its efficiency. Conclusion and outlook are 
provided in Section 8. 

2. Modelling. We consider the original model introduced in [26], which is a 
conserved gradient flow of a Swift-Hohenberg energy and serves as a model system for 
a regular periodic wave-like order-parameter field that can be interpreted as particle 
density. The Swift-Hohenberg energy is given here in a simplified form: 

FW = jf ^ + ^(r + (l+A) a )*da, (2.1) 

where the order-parameter field ip describes the deviation from a reference density, the 
parameter r can be related to the temperature of the system and fl C R m , m— 1, 2, 3 
is the spatial domain. According to the notation in [65] we consider the 7J _1 -gradient 
flow of F, the PFC2-model: 


dtip = A 


6F[ip] 

Sip 


( 2 . 2 ) 


respective a Wasserstein gradient flow [43] of F, the PFCl-model, as a generalization 
of (2.2): 


d t iP = V- (V»+V^M), (2.3) 

with a mobility coefficient ip + = ip — ip m i n > 0 with the lower bound 1 ip m i n = —1-5. 
By calculus of variations and splitting of higher order derivatives, we can find a set 
of second order equations, which will be analyzed in this paper: 

= ip 3 + (1 + r)ip + 2A ip + Aw in H x [0, T] 
d t ip = V • (i/> + V/x) (2.4) 

ui = Aip, 


for a time interval [0,T] and subject to initial condition ip(t = 0) = ipo in II and 
boundary conditions on dfl, e.g. homogeneous Neumann boundary conditions 

d„ip = d n tjJ = ip + d n ii = 0 on dfl. 


1 Tlie lower bound y tn iri = —1.5 is due to the scaling and shifting of the order-parameter from a 
physical density with lower bound 0. 
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3. Discrete equations. To transform the partial differential equation (2.4) into 
a system of linear equations, we discretize in space using finite elements and in time 
using a backward Euler discretization, respective a Rosenbrock discretization scheme. 

Let fl C K m be a regular domain (m = 1,2,3) with a conforming triangulation 
Th(Cl) with h = max^g^ (hx) a discretization parameter describing the maximal 
element size in the triangulation. We consider simplicial meshes, i.e. made of line 
segments in ID, triangles in 2D and tetrahedra in 3D. Let 

V h := {v G H\n) ; v\ T G P p (T), VT G T h (Sl)} 

be the corresponding finite element space, with P p (T) the space of local polynomials 
of degree < p , where we have chosen p = 1,2 in our simulations. The problem (2.4) 
in discrete weak form can be stated as: 

Find (J, h ,i/) h ,(jJ h G L 2 (0,T; 14) with ip h (t = 0) = i/j 0 G L 2 (£l), s.t. 

(p h ~ tph ~ (1 + r)iph,fih)si + (2 V'lph + Vwh, Vi4)n 

+ (dt^ h ,K)n + (^V^,V<) n (3.1) 

+ K, <) n + (V^ fe , V<)n = 0 W h , 0' h , < G 14, 

with (u, u)n := f n u ■ v dx. 

In the following let 0 = to < ti < • • • < tjv = T be a discretization of the time 
interval [0,T]. Let r ). := tk+i — tk be the timestep width in the fc-th iteration and 
ij)k = respective pk = Hh{tk) and u>k = Uh{tk) the discrete functions at time 

tk- Applying a semi-implicit Euler discretization to (3.1) results in a time and space 
discrete system of equations: 

Let V’o G L 2 (£l) be given. For k = 0,1,..., IV - 1 find /J-k+i, i’k+i, ^k+i G 14, s.t. 
a (e) ((/Xfe+i,'0fc + i,w fc+ i), (’d h ,'&h,K)) := 

{p-k+i - (1 + r)ip k +i,&h)n + (2V4fc+i + Vwfe+i, Vi?fc)n 

+ (ipk+i,^’h) n + (rfcV’fc V/Tfe+i, V??^)n 

+ (u> k+1 ,K) n + (V4 fc+ i,V<)n 

= + tyfcA) n = : <^ (e) , AAA)> e 14- 

(3.2) 

Instead of taking explicitly it is pointed out in [12], that a linearization of this 
non-linear term stabilizes the system and allows for larger timestep widths. Therefore, 
we replace (^,i9/,)n by (3i/jl4> k +i - 24|,t4)n. Thus (3.2) reads 

a((p k +i,i>k+i, Wfc+i), {$h, K, K)) '■= 

(fJ-k+i ~ (1 + r)i/jk+ i - 3il)ltp k +i,^h)n + (2V4 fc+ i + Vw fc+ i, V?4)n 

+ (V>fc+i A ) n + ('TfeV’fc V/Xfe+1, Vd'Jn (3.3) 

+ ( Wt+ iX)n + (Vk,V^)n 

= (-2^,i? h ) n + =: <F,(^,<,<)) W h ,K,< e Vi,, 

Let {(^} be a basis of 14, than we can define the system matrix A and the 
right-hand side vector b, for the linear system Ax = b, as 


A = 

A 00 

A 10 

A 01 

A 11 

A 02 ' 

A 12 

, b = 



A 20 

A 21 

A 22 
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with each block defined via 


[A t3 ] k i = a(ej<pi,enp k ), [b l ]j = (F,enpj), 


where e ? ; is the *-th Cartesian unit vector. 

Introducing the short cuts M := ((tpj, <fii)n)^ and K := ((V^-,Vy>»)n) for 
mass- and stiffness-matrix, K+(^) := . for the mobility matrix 

and for the non-linear term the short cut N(^) := ((—3we can write 
A as 


A = 


M 

rfcK+^fe) 

0 


— (1 + r)M + N(ipk) + 2K K 

M 0 

K M 


(3.4) 


We also find that b° = ((—2ipu,ipj)n) b 1 = M if’ and b 2 = 0. Using this, we can 
define a new matrix B := KM _1 K to decouple the first two equations from the last 
equation, i.e. 


A ,_ [ M -(l + r)M + N(^ fc )+2K-B] 

A - [r k K + (^ k ) M J ' 

With x = it± k+v ± k+v u k+1 ) T ,x' = (M fc+1 ,^ +1 ) T ,b' = (b°,b 1 ) T , where the dis¬ 
crete coefficient vectors correspond to a discretization with the same basis functions 
as the matrices, i.e. 


iph = E i/>(i)<Pi with coefficients i/j = (^(i))*; 

i 


and /i, uj_ in a same manner, we have 

(3.3) » Ax = b « AV = b', Mw fc+1 = -K^ . 

The reduced system can be seen as a discretization of a partial differential equation 
including the Bi-Laplacian, i.e. 

d t ip = V • (ijj + V/z), with n = ip 3 + (1 + r)ip + 2A ip + A 2 ip. 

In the following, we will drop the underscore for the coefficient vectors for ease of 
reading. 

3.1. Rosenbrock time-discretization. To obtain a time discretization with 
high accuracy and stability with an easy step size control, we replace the discretization 
(3.2), respective (3.3), by an embedded Rosenbrock time-discretization scheme, see 
e.g. [36, 46, 55, 41, 42], 

Therefore consider the abstract general form of a differential algebraic equation 


Md f x = F[x], 

with a linear (mass-)operator M and a (non-linear) differential operator F. Using 
the notation JJ]p(x)[y] := g^F[x + ey] | e _ 0 for the Gateaux derivative of F at x in the 
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direction y, one can write down a general Rosenbrock scheme 


—My? - JF(x fc )[yf] = F[x?] + V ^My* 

Tfe7 Z —' Tk 

0 =1 

for * = 1 ,..., s 

x? = x fc + a^y^ (* th stage solution) 
j=i 

X fc+1 = X fc + £ m jy *, x fc+1 = x fc + J > 

t=i j=i 


(3.6) 


(3.7) 


with coefficients 7 , a ^. Cjj, m,. r'rii and timestep 77. The coefficients m^ and m; build 
up linear-combinations of the intermediate solutions of two different orders. This 
can be used in order to estimate the timestep error and control the timestep width. 
Details about stepsize control can be found in [36, 46]. The coefficients used for the 
PFC equation are based on the RosSPw scheme [55] and are listed in Table 3.1. This 
W-method has 3 internal stages, i.e. s = 3, and is strongly A-stable. As Rosenbrock- 
method it is of order 3. It avoids order reduction when applied to semidiscretized 
parabolic PDEs and is thus applicable to our equations. 


7 = 0.78867513459481287 

C 11 = -C 22 = 7 
c 2 i = —2.53589838486225 
c 3 i = -1.62740473580836 
c 32 = -0.274519052838329 
c 33 = -0.0528312163512967 

021 = 2 

022 = 1.57735026918963 
oar = 0.633974596215561 
033 = 0.5 

mi = 1.63397459621556 
m 2 = 0.294228634059948 
m 3 = 1.07179676972449 

mi = 1.99444650053487 
m 2 = 0.654700538379252 
m 3 = m 3 


Table 3.1 

A set of coefficients for the Ros3Pw Rosenbrock scheme translated into the modified form of 
the Rosenbrock method used in (3.6) . All coefficients not given explicitly are set to zero. 


In case of the PFC system (2.4) we have x = (/ i,,ip,tjj) T and M = diag(0,1,0). 
The functional F applied to x is given by 



—ft + (1 + r)ip + 2A ip + Aw 


ip 3 

F[x] = 

0 

+ 

V ■ (V>+V/z) 


—to + Atp 


0 


F Lin [x] 


For the Jacobian of F in the direction y = (dfi, dtp, duj) T we find 


Jp(x)[y] = F Lin [y] 4 

(assuming ip + = ip - ip min ) 


Zip 2 dip 

V • (ip+Vdft) + V • {dipd^{ip + )V ft) 

0 


= F L in[y] 


3 ip 2 dip 

V • (ip+Vdft) + V ■ (#V/x) 
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By multiplication with test functions •& = (■&, •&', ^") T and integration over f2, we can 
derive a weak form of equation (3.6): 

For i = find y£ G ( L 2 (0,T ; V h )) 3 , s.t. 

1 * _1 

— (Myf,^) n -JI F (x fc )[y?^]=F(xf)[t9] + V^(My J fc ,i9)a W> G (^) 3 , (3.9) 
Tfe7 “ T k 

with the linear form F(-)[•]: 




= (-Ai + (l + r)V»,<?)n-(2VV> + Va;,Vi?)n 
+ (^ 3 ,^)n-(V' + VM:V^ , )n 
=: F Lin (x)[0] + (^ 3 ,^)n - ty+V/i, Vtf')n 


(w,0")n 


(V^,V0")n 


and the bi-linear form J F (•)[•, ■]: 

l F (x)[y,i9] = F Lin (y)[tf] + (3ip 2 di/j,'d)n - (ip + S/dfi + d^V/x, Vi?')n. 


Using the definitions of the elementary matrices M, K, K + and N, as above and 
introducing F(/z) := ((^ V/i. Vifii)nj , we can write the Rosenbrock discretization in 
matrix form for the i-th stage iteration: 


M -(l + r)M + 2K + N(^ fe ) K 
T k K+(V>fe) iM + r fc F(/i fe ) 0 

0 KM 




= bf 


(3.10) 


with bf the assembling of the right-hand side of (3.9), with a factor Tk multiplied 
to the second component. The system matrix A R in each stage of one Rosenbrock 
time iteration is very similar to the matrix derived for the simple backward Euler 
discretization in (3.4), up to a factor f in front of a mass matrix and the derivative 
of the mobility term F. The latter can be simplified in case of the PFC2 model (2.2), 
where F = 0 and K + = K. 

4. Precondition the linear systems. To solve the linear system Ax = b, re¬ 
spective A R y = b fi , linear solvers must be applied. As direct solvers, like UMFPACK 
[21], MUMPS [6] or SuplerLU_DIST [47] suffer from fast increase of memory require¬ 
ments and bad scaling properties for massively parallel problems, iterative solution 
methods are required. The system matrix A, respective A^, is non-symmetric, non¬ 
positive definite and non-normal, which restricts the choice of applicable solvers. We 
here use a GMRES algorithm [59], respectively the flexible variant FGMRES [58], to 
allow for preconditioners with (non-linear) iterative inner solvers, like a CG method. 

Instead of solving the linear system Ax = b, we consider the modified system 
AP _1 (Px) = b, i.e. a right preconditioning of the matrix A. A natural requirement 
for the preconditioner P is, that it should be simple and fast to solve P -1 v for 
arbitrary vectors v, since it is applied to the Krylov basisvectors in each iteration of 
the (F)GMRES method. 

We propose a block-preconditioner P for the 2x2 upper left block matrix A' of 
A based on an approach similar to a preconditioner developed for the Cahn-Hilliard 
equation [17]. Therefore, we first simplify the matrix A', respective the corresponding 
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reduced system A n ' of A R , by considering a fixed timestep T& = r and using a 
constant mobility approximation, i.e. K + ss M 0 K, with M 0 = (ip + ) the mean of the 
mobility coefficient ip + , and F = 0. For simplicity, we develop the preconditioner 
for the case Mq = 1 and 7=1 only. For small timestep widths r the semi-implicit 
Euler time-discretization (3.2) is a good approximation of (3.10), so we neglect the 
non-linear term N (ip). What remains is the reduced system 


A" := 


M 

rK 


— (1 + r)M + 2K — B 

M 


By adding a small perturbation to the diagonal of A", we can find a matrix having 
an explicit triangular block-factorization. This matrix we propose as a preconditioner 
for the original matrix A': 


M 2K B 


M 0 

' I M _ 1 (2K — B) 

rK M — <5K + dB 


rK M + <5K 

0 M -1 (M — 2<5K + <5B) 


(4.1) 


with S := yfr. In each (F)GMRES iteration, the preconditioner is applied to a vector 
(bo,bi) T , that means solving the linear system Px = b, in four steps: 

(1) My 0 = b 0 (2) (M + 5K) yi = b, - rK yi 

(3) (M - 25K + <5B) Xl = M yi (4) x 0 = y 0 + j(yi - Xl ). 

Since the overall system matrix A has a third component, that was removed for the 
construction of the preconditioner, the third component b 2 of the vector has to be 
preconditioned as well. This can be performed by solving: (5) M X 2 = b 2 — Kxy. 

In step (3) we have to solve 

Sxi := (M-2(5K + (5KM” 1 K) Xl = Myi, (4.2) 


which requires special care, as forming the matrix S explicitly is no option, as the 
inverse of the mass matrix M is dense and thus the matrix S as well. In the following 
subsections we give two approximations to solve this problem. 

4.1. Diagonal approximation of the mass matrix. Approximating the mass 
matrix by a diagonal matrix leads to a sparse approximation of S. Using the ansatz 
M ^ 1 « diag(M) 1 =: M ^ 1 the matrix S can be approximated by 

S D := (M - 2<5K + SKM^K). 

By estimating the eigenvalues of the generalized eigenvalue problem AS^x = Sx we 
show, similar as in [16], that the proposed matrix is a good approximation. 

Lemma 1. The eigenvalues A of the generalized eigenvalue problem AS^x = Sx 
are bounded by bounds of the eigenvalues g of the generalized eigenvalue problem 
/jM^y = My for mass-matrix and diagonal approximation of the mass-matrix. 

Proof. We follow the argumentation of [16, Section 3.2]. 

Using the matrices D := and K := M^5KM~ 2 we can reformulate 

the eigenvalue problem ASr>x = Sx as 

AM? (I - 26K + (5KDK)M5 X = Ms (I — 2SK + dKK)M^x. (4.3) 
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Multiplying from the left with x T , dividing by || 1VC 2 x || 2 and defining the normalized 
vector y := M 5 x/||M 2 x|| results in a scalar equation for A: 

A(1 — 28k + 8k 2 d ) = 1 — 2 5k + Sk 2 , 

with the Rayleigh quotients k = y T Ky/(y T y) and d = y T Dy/(y T y). Assuming 
that (1 — 2 5k + 5k 2 d) ^Owe arrive at 

x _ 1-2 5k + 5k 2 
~ l-25k + 5k 2 d , 

where the difference in the highest order terms of the rational function is the factor 
d. From the definition of d and D, bounds are given by the bounds of the eigenvalues 
of = Mu. □ 

In [67] concrete values are provided for linear and quadratic Lagrangian finite 
elements on triangles and linear Lagrangian elements on tetrahedra. For the latter, 
the bound d £ [0.3924,2.5] translates directly to the bound for A, i.e. A £ [0.3924, 2.5], 
and thus S d provides a reasonable approximation of S. 

Remark 1. Other diagonal approximations based on lumped mass matrices could 
also be used, which however would lead to different eigenvalue bounds. 

4.2. Relation to a Cahn-Hilliard system. An alternative to the diagonal ap¬ 
proximation can be achieved by using the similarity of step (3) in the preconditioning 
with the discretization of a Cahn-Hilliard equation [18, 17]. This equation can be 
written using higher order derivatives: 

d t c = A (c 3 ) — Ac — ? 7 A 2 c 

with 77 a parameter related to the interface thickness. For an Euler discretization in 
time with timestep width t' and finite element discretization in space as above, we 
find the discrete equation 

(M - t'K + r'rjB - r'N'(c fc ))c fc+1 = Mc k . 

Setting 77 := -i and t' := 25, and neglecting the Jacobian operator N' we recover 
equation (4.2). A preconditioner for the Cahn-Hilliard equation, see [17, 16, 8 ] thus 
might help to solve the equation in step (3), which we rewrite as a block system 

' M M - 77 K 
r'K M 

with Schur complement S. Using the proposed inner preconditioner Aq of [17, p. 13]: 

' M -77K 

r'K M + 2 v / t 7 ?/kJ ’ 

with Schur complement Sqh := M + 2 v / t 7 ?/K+t , ? 7 KM~ 1 K as a direct approximation 
of equation (4.4), respective (4.2), i.e. 

S CH xi = (M + 2v / JK + JKM^ 1 K)xi = Myi (4.5) 

we arrive at a simple two step procedure for step (3): 

(3.1) (M + V8K)z = Myi (3.2) (M + v / JK)xi = Mz. 









Lemma 2. The eigenvalues A of the generalized eigenvalue problem AS ctx = Sx 
satisfy A £ [(1 — y/6)/ 2,1], 

Proof We follow the proof of [52, Theorem 4] and denote by A the eigenvalue 
of SpjljS with the corresponding eigenvector x. We have M symmetric and positive 
definite and hence I + \/(5M _1 K positive definite and thus invertible. 

SqhSx = Ax 

=> (M + 2v^K + <5KM - 1 K) -1 (M - 2<5K + 5KM _ 1 K)x = Ax 
=> (I + v^M- 1 K )" 2 (I - 2(5M _1 K + <5(M _ 1 K) 2 )x = Ax. 

Thus, for each eigenvalue p of M _1 K we have p £ R>o and 


X(p) — (p 2 + 2 Sp + 6)(p + y/~8) 2 

an eigenvalue of Sq^S and since M 1 K is similar to M 1 / 2 M 'KM x l 2 that is sym¬ 
metric, all eigenvalues are determined. 

With algebraic arguments and y/S > 0 we find 


•M p) < 


p 2 + (V6) 2 

(,h + VS ) 2 


< 1 


and VA = 0 for p, S \ 0. This leads to the lower bound < A (p). □ 

With this Lemma we can argue that Sch provides a good approximation of S at 
least for small timestep widths r = 6 2 < 1 . 

We can write the matrix P = P(S) in terms of the Schur complement matrix S: 


P(S) 


M (5 -1 (M — S) 

rK SK + S 


(4.6) 


Inserting Sch instead of S gives the precondition-matrix for the Cahn-Hilliard ap¬ 
proximation 


Pch := P(Sch) 


M -~2y/6K - B 

rK M + {6 + 2y/S)K. + 5B 


5. Convergence analysis of the Krylov-subspace method. To analyze the 
proposed preconditioners for the GMRES algorithm, we have a look at the norm of 
the residuals rfc(A) = b — Ax*, of the approximate solution x* obtained in the fc-th 
step of the GMRES algorithm. In our studies, we are interested in estimates of the 
residual norm of the form 


IHh 

11 r 0 112 


min 

pGiifc 


lb( A ) r oll 2 

lko || 2 


< mm ||p(A)[| 2 


(5.1) 


with II*, := {p £ P* : p(0) = 1} and r 0 the initial residual. The right-hand side corre¬ 
sponds to an ideal-GMRES bound that excludes the influence of the initial residual. 
In order to get an idea of the convergence behavior, we have to estimate / approximate 
the right-hand side term by values that are attainable by analysis of A. Replacing A 
by AP” 1 we hope to get an improvement in the residuals. 
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A lower bound for the right-hand side of (5.1) can be found by using the spectral 
mapping theorem p(cr(A)) = a(p( A)), as 

min max |p(A)| < min ||p(A)|| 2 , (5.2) 

pen fc A6 (t(A) pen fc 

see [64, 23] and an upper bound can be stated by finding a set 5(A) C C associated 
with A, so that 


min ||p(A)|| 2 < C min max |p(A)|, (5.3) 

pen k pen fc agS(A) 


where C is a constant that depends on the condition number of the eigenvector matrix, 
the e-pseudospectra of A respective on the fields of values of A. 

Both estimates contain the min-max value of p( A). In [56, 23] it is shown, that 
the limit 


lim 

k—>oo 


1 1 / k 

min max |p(A)| =: ps 

pGlIfc A£<S 


exists, where ps is called the estimated asymptotic convergence factor related to the 
set S. Thus, for large k we expect a behavior for the right-hand side of (5.1) like 

Pl{ A) ^ min |b(A)|| 2 < Cp k s{A) . 


The tilde indicates that this estimate only holds in the limit k —> oo. 

In the next two sections, we will summarize known results on how to obtain the 
asymptotic convergence factors ps and the constant C in the approximation of the 
relative residual bound. 


5.1. The convergence prefactor. The constant C plays an important role in 
the case of non-normal matrices, as pointed out by [30, 64], and can dominate the 
convergence in the first iterations. It is shown in Section 6 that the linear part of 
the operator matrix related to A is non-normal and also the preconditioned operator 
related to Q := AP 1 is non-normal. Thus, we have to take a look at this constant. 

An estimate of the convergence constant, applicable for general non-normal ma¬ 
trices, is related to the e-pseudospectrum <r e ( A) of the matrix A. This can be defined 
by the spectrum of a perturbed matrix [64, 30] 

ct £ (A) := {A £ C | z £ <t(A + E), ||E|| 2 5 e} ■ 

Let T e := da e be the boundary of cr e , respective an union of Jordan curves ap¬ 
proximating the boundary, then 


|rj |r c | 

min |b(A)|| 2 < min max |p(A)| < min max |p(A)|, (5.4) 

pe 11*, 27re pen k Aer, 27re p& n fc Aeo- e (A) 

and thus C = with |T e | the length of the curve T e [64]. This estimate is approx¬ 
imated, using the asymptotic convergence factor for large k , by 

min ||p(A)|| 2 < ^p \\. < ^^.(A)- (5-5) 

This constant gives a first insight into the convergence behavior of the GMRES 
method for the PFC matrix A respective the preconditioned matrix Q. 
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5.2. The asymptotic convergence factor. The asymptotic convergence fac¬ 
tor ps, where S is a set in the complex plane, e.g. S = <r(A), or S = <t € (A), can 
be estimated by means of potential theory [23, 44]. Therefore, we have to construct 
a conformal mapping $ : C —> C of the exterior of S to the exterior of the unit disk 
with <h(oo) = oo. We assume that S is connected. Otherwise, we will take a slightly 
larger connected set. Having S C C \ {0} the convergence factor is then given by 


\m\' 


(5.6) 


Let S = [a, (3] be a real interval with 0 < a < (3 and «:=■£, then a conformal 
mapping from the exterior of the interval to the exterior of the unit circle is given by 




2z — K — 1 — 2~J z 2 — (k + l)z + K 

K — 1 


(5.7) 


see [23] 2 , and gives the asymptotic convergence factor 


P[<X,1 3] 


\[k — 1 

\/n +1 


(5.8) 


that is a well known convergence bound for the CG method for symmetric positive 
definite matrices with n = the spectral condition number of the matrix A. In 

•^min 

case of non-normal matrices, the value n is not necessarily connected to the matrix 
condition number. 

In the next Section, we will apply the given estimates for the asymptotic conver¬ 
gence factor and for the constant C to the Fourier transform of the operators that 
define the PFC equation, in order to get an estimate of the behavior of the GMRES 
solver. 


6. Spectral analysis of the preconditioner. We analyze the properties and 
quality of the proposed preconditioner by means of a Fourier analysis. We therefore 
consider an unbounded, respective periodic domain VL and introduce continuous op¬ 
erators A, Ao and P associated with the linear part F^n of (3.8), of the linear part of 
the 6th order non-splitted version of (2.2) and the preconditioner, for xjj + = 1: 


p — (1 + r)ip — 2Aip — A ui 
— tA/z + if} 
u> — Ai/> 

Ao[0] := ijj — rA((l + r)ip + 2Aip + A 2 ip ), 
with x = The operator, that represents the preconditioner reads 

p — 2A'i j) — A 2 ip 
—rAp if) — 8Axj) + 8A 2 xp 
ui — A xjj 


Using the representation of P(S) in (4.6), we can also formulate the operator that 
determines the Cahn-Hilliard approximation of P by inserting SciP 


P C h[x] := 


p + 2\ZSAx/j — A 2 xp 
-tA/i + xj} + (8 - f 2\/S)Axp + 5A 2 xj) 
u> — Axj) 


2 In [23] the sign of the square-root is wrong and thus, the exterior of the interval is mapped to 
the interior of the unit circle. In formula (5.7) this has been corrected. 


11 














We denote by k = (fc 1; fc 2 , k 3 ) the wave vector with k 2 = fc 2 + fcf+ fc 2 . The Fourier 
transform of a function u = u(r) will be denoted by u = u(k) and is defined as 

T : u(r) i->- u( k) = f e _l(k ' r ^u(r) dr. 

Jr 3 

Using the inverse Fourier transform, the operators A,Ao,P and Pch applied to x, 
respective ip, can be expressed as A[x] = F~ l {Aft), A 0 [^] = J 7 ~ 1 (Aoip), P[x] = 
J r ~ 1 (PSt) and Pch[x] = .F _ 1 (Pchx), with x = (p,,ip,Q) and Aq,A,V and Pch the 
symbols of Ao, A, P and Pch, respectively. These symbols can be written in terms of 
the wave vector k: 


Ax = 


1 

rk 2 

0 


-(1 + r) + 2 k 2 

1 


k 2 



Aoip = 

Px = 


Pchx = 


[l+T 

[(1 + r)k 2 - 2k 4 

+ 

k 6 ])^, 

' 1 

2k 2 - k 4 

O' 

(P\ 

rk 2 

1 - <5k 2 + 5k 4 

0 


_ 0 

k 2 

1 

W 

' 1 

-2v^k 2 - 

k 4 

o' 


rk 2 1 + (6 + 2VS)k 2 + 5k 4 0 

0 k 2 1 



( 6 . 1 ) 

( 6 . 2 ) 

(6.3) 

(6.4) 


In Figure 6.1 the eigenvalue symbol curves of A restricted to a bounded range of 
frequencies, together with the distribution of eigenvalues of an assembled matrix 3 A, 
using quadratic finite elements on a periodic tetrahedral mesh with grid size h = 7 t/4 , 
is shown. The qualitative distribution of the eigenvalues is similar for symbol curves 
and assembled matrices, and changes as the timestep width increases. 

For small r, the origin is excluded by the Y-shape profile of the spectrum. In¬ 
creasing r leads to a surrounding of the origin. This does not necessarily imply a bad 
convergence behavior. 

6.1. Critical timestep width. For larger timestep widths r we can even get the 
eigenvalue zero in the continuous spectrum, i.e. the time-discretization gets unstable. 
In the following theorem this is analyzed in detail and a modification is proposed, 
that shifts this critical timestep width limit. This modification will be used in the 
rest of the paper 

Lemma 3. Let Ao be given as in (6.2) and r < 0 then the spectrum <r(Ao) 
contains zero in case of the critical timestep width 


T > T* 


27 

2{y/a — l)(y/a + 2) 2 ’ 


(6.5) 


with a = 1 — 3 r. 


3 Since a finite element mass-matrix has eigenvalues far from the continuous eigenvalue 1, de¬ 
pending on the finite elements used, the grid size and the connectivity of the mesh, the overall 
spectrum of A (that contains mass-matrices on the diagonal) is shifted on the real axis. In order 
to compare the continuous and the discrete spectrum, we have therefore considered the diagonal 
preconditioned matrix A = diag(A) -X A that is a blockwise diagonal scaling by the inverse of the 
diagonal of mass-matrices. 
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o eigenvalues of an assembled FEM matrix 
- continuous operator spectrum 



Fig. 6.1. Eigenvalues of the diagonally preconditioned finite element matrix A = diag(A) —1 A ? 
i.e. a discretization of the continuous operator A multiplied with the inverse of its diagonal, and the 
three eigenvalues of the symbol A visualized as restricted symbol curves. Left: spectrum for timestep 
width r = 0.1, Right: spectrum for timestep width r = 1 


Let The spectrum of the modified operator Aq, given by 

A 0 := (1 + t[( 3if 2 + 1 + r)k 2 - 2k 4 + k 6 ]), (6.6) 

contains zero only in the case of ip' 2 < A/. Then the critical timestep width is given 
by (6.5) with a = 1 — 3r — 9ip 2 . 

Remark 2. The modified operator An can be derived by linearizing ip 3 around a 
constant reference density ip: 

ip 3 « 3ip 2 ip — 2 ip 3 . 

Adding this as an approximation of the non-linear term to the system (6.2) leads to 
the operator symbol (6.6). 

Remark 3. If we take ip as the constant mean value of ip over fl, with r has the 
physical meaning of an undercooling of the system, then the relation \ip\ = \J —r /3 
is related to the solid-liquid transition in the phase-diagram of the PFC model, i.e. 
\ip\ > \J—r/3 leads to stable constant solutions, interpreted as liquid phase, and \ip\ < 
sf—r /3 leads to an instability of the constant phase, interpreted as crystalline state. 
An analysis of the stability condition can be found in [20, 27] upon other. 

Remark 4. In [68, 38] an unconditionally stable discretization is provided that 
changes the structure of the matrix, i.e. the negative 2k 4 term is moved to the right- 
hand side of the equation. In order to analyze also the Rosenbrock scheme we can not 
take the same modification into account. The modification shown here is a bit similar 
to the stabilization proposed in [29], but the authors have added a higher order term 
Ck 4 instead of the lower order term C'k 2 in (6.6). 

Proof. (Lemma 3). We analyze the eigenvalues of An and get the eigenvalues of 
Aq as a special case for ip = 0. The eigenvalue symbol Aq gets zero whenever 


(3ip 2 + 1 + r)k 2 — 2k 4 + k 6 

The minimal r 6 R>ch denoted by r*, that fulfils this equality is reached at 

.9 2 1 /-- 2 1 _ 

k~ =- +-fl-3r-9ip 2 =: - + -s/a. 


13 
















Inserting this into r gives 


r 


* 



We have a > 0 for |i/;| > ^y/l — 3 r and r* > 0 4=> a > 1 4=> ^ by simple 

algebraic calculations. 0 

On account of this zero eigenvalue, we restrict the spectral analysis to small 
timestep widths r. For r = —0.35, as in the numerical examples below, we get the 
timestep width bound 0 < r < r* ~ 2.6548 for the operator Aq and with ?/> = —0.34 
we have 0 < r < r* « 312.25 for Aq, hence a much larger upper bound. In the 
following, we will use the modified symbols for all further calculations, i.e. 


A = 


1 

rk 2 

0 


— (3 ip 2 + 1 + r) + 2k 2 k 2 
1 0 

k 2 1 


(6.7) 


and remove the hat symbol for simplicity, i.e. A —> A, Aq —> Aq. 

Calculating the eigenvalues of the preconditioner symbol Q := AT > ~ 1 , respective 
Qch := directly gives the sets 


c(Q) 

ct(Qch) 



r(k 6 - 2k 4 + (3V> 2 + 1 + r)k 2 ) + 1 
rk 6 + (i Jt — 2r)k 4 — s/rk 2 + 1 

r(k 6 — 2k 4 + (3t/> 2 + 1 + r)k 2 ) + 1 
rk 6 + (i/t + 2r 3 / 4 )k 4 + (a/t + 2r 1 / 4 )k 2 + 1 


k G 


( 6 . 8 ) 

(6.9) 


with values all on the real axis (for r > 0). Similar to the analysis of Aq we get 
a critical timestep width, i.e. eigenvalues zero, for r > r*. The denominator of the 
third eigenvalue of ct(Qch) is strictly positive, but the denominator of cr(<2) can reach 
zero. This would lead to bad convergence behavior, since divergence of this eigenvalue 
would lead to divergence of the asymptotic convergence factor in (5.8). 

The critical timestep width, denoted by , that allows a denominator with value 
zero is given by r = (—k 4 + 2k 2 )^ 2 that is minimal positive for |k| = 1 and gives 
A = 1. Thus, for the preconditioner V we have to restrict the timestep width to 
r G (OjT 11 ) C (0,t*). This restriction is not necessary for the preconditioner Vcn- 

6.2. The asymptotic convergence factor. Since the third eigenvalue of Q 
in (6.8), respective Qch in (6.9), is a real interval C M+, for r in the feasible range 
(0, min(r l! ,T*)), we can use formula (5.8) to estimate an asymptotic convergence factor 
for the lower bound on min pe n fc ||p(Q*)||- For fixed r = —0.35 and various values of Tp 
minimum and maximum of (6.8) and (6.9) are calculated numerically. Formula (5.8) 
thus gives the corresponding estimated asymptotic convergence factor (see left plot 
of Figure 6.2). For step widths r less than 1 we have the lowest convergence factor 
for the operator Q and the largest for the original operator Aq. The operator Qch 
has slightly greater convergence factor than Q but is more stable with respect to an 
increase in timestep width. 

The stabilization term 3^ 2 added to Aq and A in (6.6), respective (6.7), influences 
the convergence factor of Aq and Q only slightly, but the convergence factor of Qch 
is improved a lot, i.e. the critical timestep width is shifted towards positive infinity. 

The upper bounds on the convergence factor are analyzed in the next Section. 
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Fig. 6.2. Left: Asymptotic convergence factor for operators Ao,Q and Qch■ In dashed lines, 
the dependence on the mean density ip-modification (6.6) is shown. Right: Comparison of the 
convergence factor related to the spectrum and e-pseudospectrum is shown. This corresponds to 
lower and upper bounds of the actual asymptotic convergence factor. 


6.3. Analysis of the pseudospectrum. As it can be seen by simple calcula¬ 
tions, the symbol A is non-normal: 

(A T A - AA t ) 2 0 = k 2 ((3<p + 1 + r) - 2k 2 ) ^ 0 

for a matrix entry at row 2 and column 0. For slightly more complex calculations it 
can be shown, that also Q := AV -1 and Qch := AP^h are non-normal: 

(Q t Q - QQ T ) 2j2 = (QchQch - QchQ5h) 2 , 2 = k 4 /0. 

For non-normal matrices we have to analyze the e-pseudospectrum in order to get 
an estimate of convergence bounds for the GMRES method, as pointed out in Section 
5.1. 

Using the Matlab Toolbox Eigtool provided by [69] we can calculate the pseu¬ 
dospectra <7 e and approximations r e of its boundaries with single closed Jordan curves 
for all wave-numbers ki 6 [0,fc max ]. The maximal frequency used in the calculations 
is related to the grid size h of the corresponding triangulation, fc max = In all the 
numerical examples below, we have used a grid size h = j that can resolve all the 
structures sufficiently well, thus we get fc max = 4 that leads to |k| max = y/mk max . 

The e-pseudospectrum of Q and Qch can be seen in Figure 6.3 for various values 
of e. The pseudospectrum of Qch gets closer to the origin than that of Q, since the 
eigenvalues get closer to the origin as well. The overall structure of the pseudospectra 
is very similar. 

For the convergence factor corresponding to the pseudospectra, we compute the 
inverse conformal map \F = <F _1 of the exterior of the unit disk to the exterior of 
a polygon So approximating the set 5 = cr e , by the Schwarz-Christoffel formula, 
using the SC Matlab Toolbox [25, 24], A visualization of the inverse map 'F for the 
pseudospectrum of A and Q can be found in Figure 6.4. 

Evaluating the asymptotic convergence factor depending on the e-pseudospectrunr 
of the matrices is visualized in Figure 6.5. The calculation is performed for fixed 
timestep width r = 0.1 and parameters r = —0.35 and ip = 0. 
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Fig. 6.3. e-Pseudospectra of the preconditioned matrices for various values of e. Left: cr e (Q), 
Right: cr e (Q,CH )> for r = 0.1 and k in the restricted range [0, |k| m ax]. The dashed lines correspond 
to the zero-axis and indicate the origin. 




Fig. 6.4. Inverse conformal map 4/ of the unit disk to the exterior of a polygon enclosing the 
e-pseudospectrum of A, respective Q for one e and a restricted range of frequencies. For both plots 
the e value is chosen small enough to have 0 in the exterior of the peudospectrum. 


Increasing e increases the radius of the sphere like shape around the point 1. For a 
simple disc the convergence factor is proportional to the radius (see [23]), thus we find 
increasing convergence factors also for our disc with the tooth. When e gets too large 
the pseudospectrum may contain the origin that would lead to useless convergence 
bounds, since then > 1 in (5.6). If e gets too small, the convergence constant C in 
(5.4) is growing rapidly, since |r e | is bounded from below by the eigenvalue interval 
length, i.e. |r e | < 2(/3 — a) = 2((max(cr(Q)), min(cr(Q))), and we divide by e. Thus, 
the estimates also are not meaningful in the limit e —> 0. 

An evaluation of the constant C for various values e can be found in the left plot 
of Figure 6.5. It is a log-log plot with constants in the range [10 2 ,10 4 ]. 

For all e > 0 the upper bound (5.4) is valid, so we have chosen e = 10 -3 and 
plotted the resulting estimated asymptotic convergence factor in relation to the lower 
bound, i.e. the convergence factor corresponding to the spectrum of the matrices, in 
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the right plot of Figure 6.2. The upper bound is just slightly above the lower bound. 
Thus, we have a convergence factor for Q that is in the range 0.2 — 0.3 (for r = 0.1 
and r = —0.35) and for the matrix Qch in the range 0.45 — 0.55, that is much lower 
than the lower bound of the convergence factor of Ao (approximately 0.99). 

So both, the preconditioner V and 'Pch improves the asymptotic convergence 
factor a lot and we expect fast convergence also in the case of discretized matrices. 




Fig. 6.5. Left: Estimated convergence prefactor C := for the matrix Q = A‘P 1 and 
Qch — •A'Pqa with r = —0.35 and r = 0.12, plotted in logarithmic scale for e and C. Right: 
Estimated asymptotic convergence factor Pa e (Q) — (Q) (0)| 1 an d fPJQcn) analogously. The 

legend in the left plot is valid also for the right plot. 


7. Numerical studies. We now demonstrate the properties of the precondi¬ 
tioner numerically. We consider a simple crystallization problem in 2D and 3D, start¬ 
ing with an initial grain in a corner of a rectangular domain. The solution of the 
PFC equation in the crystalline phase is a periodic wave-like field with specific wave 
length and amplitude. In [27, 40] a single mode approximation for the PFC equation 
in 2D and 3D is provided. These approximations show a wave length of d := 47r /y/3, 
corresponding to a lattice spacing in a hexagonal crystal in 2D and a body-centered 
cubic (BCC) crystal in 3D. We define the domain H as a rectangle/cube with edge 
length, a multiple of the lattice spacing: D = [N ■ d] 2 ’ 3 , with N £ N>o- Discretizing 
one wave with 10 gridpoints leads to a sufficient resolution. Our grid size therefore is 
h = ss j throughout the numerical calculations. We use regular simplicial meshes, 
with h corresponding to the length of an edge of a simplex for linear elements and 
twice its length for quadratic elements to guarantee the same number of degrees of 
freedom (DOFs) within one wave. 

7.1. General problem setting and results. As system parameters we have 
chosen values corresponding to a coexistence of liquid and crystalline phases: 2D 
(tp = —0.35, r = —0.35), 3D (tp = —0.34, r = —0.3). Both parameter sets are 
stable for large timestep widths, with respect to Lemma 3. In Figure 7.1 snapshots of 
the coexistence regime of liquid and crystal is shown for two and three dimensional 
calculations. 

In the steps (1), (2), (3) and (5) of the preconditioner solution procedure, lin¬ 
ear systems have to be solved. For this task we have chosen iterative solvers with 
standard preconditioners and parameters as listed in Table 7.1. The PFC equation is 
implemented in the finite element framework AMDiS [66, 57] using the linear algebra 
backend MTL4 [33, 22, 34] in sequential calculations and PETSc [15] for parallel cal¬ 
culation for the block-preconditioner (4.1) and the inner iterative solvers. As outer 
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Fig. 7.1. Intermediate state of growing crystal, starting from one corner of the domain. Shown 
is the order parameter field if. Left: D = [20 d] 2 , number of DOFs: 263,169, calculated on 1 
processor; Right: D = [12c£] 3 , number of DOFs: 101,255,4%7, calculated on 3.456 processors. 


solver, a FGMRES method is used with restart parameter 30 and modified Gram- 
Schmidt orthogonalization procedure. The spatial discretization is done using La¬ 
grange elements of polynomial degree p — 1,2 and as time discretization the implicit 
Euler or the described Rosenbrock scheme is used. 


precon. steps 

matrix 

solver 

precond. 

rel. tolerance 

(1).(5) 

M 

PCG 

diag 

10" 3 

(2) 

M + <5K 

PCG 

diag 

1CT 3 

(3.1), (3.2) 

M +v^K 

PCG 

diag 

1CT 3 

(3) 

M - 2SK + (5KM^K 

PCG 

diag 

20 (iter.) 


Table 7.1 


Parameters for the inner solvers of the preconditioner with Cahn-Hilliard approximation S ch 
and in the last line for the diagonal approximation Sp of the matrix S. ‘PCG’ is the shortcut 
for preconditioned conjugate gradient method. The preconditioner named ‘diag’ indicates a Jacobi 
preconditioner. We have solved each inner system up to a relative solver tolerance given in the last 
column of the table. Only in the case of the matrix Sp it is more efficient to use a fixed number of 
iteration. 


The first numerical test compares a PFC system solved without a preconditioner 
to a system solved with the developed preconditioner. In Figure 7.2 the relative 
residual in the first timestep of a small 2D system is visualized. For increasing timestep 
widths the FGMRES solver without preconditioner (dashed lines) shows a dramatic 
increase of the number of iterations up to a nearly stagnating curve for timestep 
widths greater than 0.5. On the other hand, we see in solid lines the preconditioned 
solution procedure that is much less influenced by the timestep widths and reaches 
the final residual within 20-30 iterations. A detailed study of the influence of the 
timestep width can be found below. For larger systems, respective systems in 3D, we 
get nearly no convergence for the non-preconditioned iterations. 

Next, we consider the solution procedure of the sub-problems in detail. Table 7.2 
shows a comparison between an iterative preconditioned conjugate gradient method 
(PCG) and a direct solver, where the factorization is calculated once for each sub¬ 
problem matrix per timestep. The number of outer iterations increases, when we use 
iterative inner solvers, but the overall solution time decreases since a few PCG steps 
are faster than the application of a LU-factorization to the Krylov vectors. This holds 
true in 2D and 3D for polynomial degree 1 and 2 of the Lagrange basis functions. 

We now compare the two proposed preconditioners regarding the same problem. 
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Fig. 7.2. Relative residual of the solver iterations. Solid lines show preconditioned solver 
iterations and dashed lines iterations without a preconditioner. The systems is $2 = [d] 2 , h = ^. 




direct 

iterative 

dim. m 

poly, degree p 

time [sec] 

^iterations 

time [sec] 

^iterations 

2D 

1 

6.11 

14 

1.65 

14 

2 

5.89 

14 

3.04 

15 

3D 

1 

41.72 

17 

3.36 

18 

2 

35.24 

17 

9.62 

19 


Table 7.2 

Comparison of the number of iterations and time to solve the linear system averaged over 20 
timesteps for the preconditioner matrix S ch with timestep width r = 0.1. Sub-problems of the 
preconditioner are solved with iterative solvers as in Table 7.1 or with the direct solver UMFPACK. 
The benchmark configuration is a problem with approximately 66,000 DOFs and grid size h = 7r/4. 


Table 7.3 shows a comparison of the approximation of the sub-problem (3) by either 
the diagonal mass-matrix approximation Sp or the Cahn-Hilliard preconditioner ap¬ 
proximation Sch- In all cases, the number of outer solver iterations needed to reach 
the relative tolerance and also the time for one outer iteration is lower for the Sch 
approximation than for the approximation. 




S 

CH 



dim. m 

poly, degree p 

time [sec] 

^iterations 

time [sec] 

^iterations 

2D 

1 

1.65 

14 

2.72 

16 

2 

3.04 

15 

7.03 

20 

3D 

1 

3.36 

18 

8.14 

21 

2 

9.62 

19 

81.49 

55 


Table 7.3 

Comparison of the number of iterations and time to solve the linear system averaged over 20 
timesteps for the preconditioner with diagonal approximation Sp of S ; respective Cahn-Hilliard 
approximation S ch- Sub-problems of the preconditioner are solved with iterative solvers as in Table 
7.1. The benchmark configuration is a problem with approximately 66,000 DOFs, timestep width 
r = 0.1 and grid size of h = 7r/4. 
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In the following, we thus use the preconditioned solution method with S^h and 
PCG for the sub-problems. We analyse the dependence on the timestep width in detail 
and compare it with the theoretical predictions and show parallel scaling properties. 

7.2. Influence of timestep width. In Table 7.4 the time to solve the lin¬ 
ear system averaged over 20 timesteps is listed for various timestep widths r. All 
simulations are started from the same initial condition that is far from the stationary 
solution. It can be found that the solution time increases and also the number of outer 
solver iterations increases. In Figure 7.4, this increase in solution time is visualized 
for various parameter sets for polynomial degree and space dimension. The behaviour 
corresponds to the increase in the asymptotic convergence factor (see Figure 6.2) for 
increasing timestep widths. 



2D 

3D 

timestep width r 

time [sec] 

^iterations 

time [sec] 

^iterations 

0.01 

2.50 

13 

8.01 

17 

0.1 

3.05 

15 

9.62 

19 

1.0 

4.53 

19 

14.29 

24 

10.0 

10.81 

47 

34.94 

58 


Table 7.4 

Comparison of time to solve the linear system averaged over 20 timesteps for various timestep 
widths r for a 2D and a 3D system. The benchmark configuration is a problem with polynomial 
degree p = 2 with approximately 66,000 DOFs. 


We have analyzed whether a critical timestep width occurs in the two approxi¬ 
mations of S (see Figure 7.3). The diagonal approximation Sd is spectrally similar 
to the original preconditioner S that has shown the critical timestep width = 1. 
In the numerical calculations, Sd shows a critical value around the analytical value, 
but it varies depending on the finite element approximation of the operators. For 
linear Lagrange elements we see Ri 2 and for quadratic element r ^ ss 0.6. The 
Cahn-Hilliard approximation Sqh does not show a timestep width, where the number 
of outer iterations explodes, at least in the analyzed interval r £ [10” 3 ,10 1 ]. The 
difference in the finite element approximations is also not so pronounced as in the 
case of Sd- 

While in all previous simulations an implicit Euler discretization was used, we now 
will demonstrate the benefit of the described Rosenbrock scheme, for which the same 
preconditioner is used. Adaptive time stepping becomes of relevance, especially close 
to the stationary solution, where the timestep width needs to be increased rapidly 
to reduce the energy further. In order to allow for large timestep widths the 

iterative solver, respective preconditioner, must be stable with respect to an increase 
in this parameter. 

In Figure 7.5 the system setup and evolution for the Rosenbrock benchmark is 
shown. We use 10 initial grains randomly distributed and oriented in the domain and 
let the grains grow until a stable configuration emerges. When growing grains touch 
each other, they build grain boundaries with crystalline defects. The orientation of 
the final grain configuration is shown in the right plot of Figure 7.5 with a color coding 
with respect to an angle of the crystal cells relative to a reference orientation. 

The time evolution of the timestep width obtained by an adaptive step size con¬ 
trol and the evolution of the corresponding solver iterations is shown in the left plot 
of Figure 7.6. Small grains grow until the whole domain is covered by particles. 
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Fig. 7.3. Increase in number of outer iterations when the timestep width increases. The 
red curves (curcular dots) correspond to the diagonal approximation So of S and the blue curves 
(asterisk dots) to the Cahn-Hilliard approximation Sch- All simulations are performed in 3D in a 
domain with grid size h = 7r/4. The solid lines correspond to simulations with polynomial degree 
p = 1 and the dashed lines with polynomial degree p = 2. 


—*— 2D domain 

- p = 1 
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Fig. 7.4. Left: Solution time per timestep iteration for various timestep widths r averaged 
over 20 timesteps. Right: Number of outer iterations per timestep for various timestep widths r 
averaged over 20 timesteps. The four curves show the dependence on dimension (2D or 3D) and on 
polynomial degree p of the Lagrange basis functions. 


This happens in the time interval [0,200], where small timestep widths are required. 
From this time, the timestep width is increased a lot by the step size control since 
the solution is in a nearly stable state. The number of outer solver iterations in¬ 
creases with increasing timestep width, as expected. Timestep widths up to 18 in the 
time evolution are selected by the step size control and work fine with the proposed 
preconditioner. 

In the right plot of Figure 7.6, the relation of the obtained timestep widths to the 
solution time is given. Increasing the timestep widths increases also the solution time, 
but the increasing factor is much lower than that of the increase in timestep width, 
i.e. the slope of the curve is much lower than 1. Thus, it is advantageous to increase 
the timestep widths as much as possible to obtain an overall fast solution time. 
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Fig. 7.5. Grain growth simulation. Left: Initial grains that do not touch each other. Center: 
Grown grains with different orientation and grain boundaries. Right: Coloring of the different 
crystal orientations. The coloring fails at the boundary of the domain. 



Fig. 7.6. Left: Time series of timestep width (solid line) and outer solver iterations (dashed 
line) for a simulation using a Rosenbrock scheme with automatic step size selection. Right: Evolu¬ 
tion of the solution time for increasing timestep widths. The time is measured relative to the time 
for the minimal timestep width. The data in extracted from the simulation of the grain growth, see 
Figure 7.5. 


7.3. Parallel calculations. We now demonstrate parallel scaling properties. 
Figure 7.7 shows strong and weak scaling results. All simulations are done in 3D and 
show results for the time to solve the linear system in comparison with a minimal 
number of processors that have the same communication and memory access environ¬ 
ment. The efficiency of this strong scaling benchmark is about 0.8 0.9 depending 

on the workload per processing unit. The efficiency of the weak scaling is about 0.9 
0.95, thus slightly better than the strong scaling. 


^processors p 

total DOFs 

time [sec] 

^iterations 

48 

1,245,456 

13.62 

24 

96 

2,477,280 

13.57 

24 

192 

4,984,512 

13.83 

25 

384 

9,976,704 

14.97 

25 


Table 7.5 

Average number of iterations and solution time for weak scaling computations. 


In Table 7.5 the number of outer solver iterations for various system sizes is given. 
The calculations are performed in parallel on a mesh with constant grid size but with 
variable domain size. By increasing the number of processors respective the number 
of degrees of freedom in the system the number of solver iterations remains almost 
constant. Also the solution time changes only slightly. 
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Larger systems on up to 3.456 processors also show that the preconditioner does 
not perturb the scaling behavior of the iterative solvers. All parallel computations 
have been done on JUROPA at JSC. 
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-o- 25,900 DOFs / core 

- - 912,673 DOFs 


19,500 DOFs / core 

. ideal speedup 

- 450 r 

-ideal speedup 
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Fig. 7.7. Speedup of parallel simulation. Left: strong scaling for fixed overall number of DOFs. 
Right: weak scaling for fixed number of DOFs per processor. The parameter q is a reference number 
of processors, that is q = 24/32 for the calculation in the left plot and q = 16/48 for the right plot. 
The first number corresponds to the dashed lines and the second to the solid lines. 


7.4. Non-regular domains. While all previous benchmark problems could 
have also been simulated using spectral or finite difference methods, we now demon¬ 
strate two examples, where this is no longer possible and the advantage of the proposed 
solution approach is shown. We use parametric finite elements to solve the surface 
PFC equation [14] on a manifold. The first example considers an elastic instability of 
a growing crystal on a sphere, similar to the experimental results for colloidal crystals 
in [49]. The observed branching of the crystal minimizes the curvature induced elastic 
energy, see Figure 7.8. The second example shows a crystalline layer on a minimal 
surface, the ‘Schwarz P surface’, Figure 7.9, which might be an approach to stabilize 
such surfaces by colloidal particles, see [39]. 




Fig. 7.8. Crystalization on a sphere <Si2o(0). Left: Visualization using OVITO [61], indicating 
each wave as a colloidal particle, Right: order parameter field if. Number of DOFs: 397,584, 
calculated on 8 processors. 


8. Conclusion and outlook. In this paper we have developed a block-precon¬ 
ditioner for the Phase-Field Crystal equation. It leads to a precondition procedure 
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Fig. 7.9. Crystal structure on a ‘Schwarz P surface’. Left and Center: Visualization using 
OVITO [61] in two different perspectives, indicating each wave as a colloidal particle, Right: order 
parameter field if. Number of DOFs: 250,000, calculated on 1 processor. 


in 5 steps that can be implemented by composition of simple iterative solvers. Ad¬ 
ditionally, we have analyzed the preconditioner in Fourier-space and in numerical 
experiments. We have found a critical timestep width for the original preconditioner 
and have proposed a variant with an inner Cahn-Hilliard preconditioner, that does 
not show this timestep limit. Since most of the calculations are performed in parallel, 
a scaling study is provided, that shows, that there is no negative influence of the 
preconditioner on the scaling properties. Thus, large scale calculations in 2D and 3D 
can be performed. 

Recently extensions of the classical PFC model are published, towards liquid 
crystalline phases [54], flowing crystals [50] and more. Analyzing the preconditioner 
for these systems, that are extended by additional coupling terms, is a planed task. 

Higher order models, to describe quasicrystalline states, respective polycrystalline 
states, based on a conserved Lifshitz-Petrich model [48], respective a multimode PFC 
model [51], are introduced and lead to even worse convergence behavior in finite 
element calculations than the classical PFC model. This leads to the question of an 
effective preconditioner for these models. The basic ideas, introduced here, might be 
applicable for the corresponding discretized equations as well. 
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